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Abstract 

We  present  a  Chebyshev  multidomain  method  that  can  solve  systems  of 
hyperbolic  equations  in  conservation  form  on  an  unrestricted  quadrilateral 
subdivision  of  a  domain.  Within  each  subdomain  the  solutions  and  fluxes 
are  approximated  by  a  staggered-grid  Chebyshev  method.  Thus,  the  method 
is  unstructured  in  terms  of  the  subdomain  decomposition,  but  strongly 
structured  within  subdomains.  Communication  between  subdomains  is 
done  by  a  mortar  method  in  such  a  way  that  the  method  is  globally 
conservative.  The  method  is  applied  to  both  linear  and  non-Unear  test 
problems  and  spectral  accuracy  is  demonstrated. 
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1 .  Introduction 


In  the  first  paper  [1],  we  introduced  a  staggered  grid  Chebyshev  multidomain 
method  for  the  solution  of  inviscid  compressible  flow  problems.  TTie  grid  used  was 
analogous  to  the  fully  staggered  grids  used  in  some  finite  difference  computations  of 
compressible  flows,  e.g.  [2].  For  the  staggered  grid  spectral  method,  the  unknowns  are 
approximated  by  global  polynomials  of  degree  A^-1  in  each  space  dimension,  which  pass 
through  values  defined  at  the  Chebyshev-Gauss  quadrature  points.  The  fluxes  are 
approximated  by  polynomials  of  degree  N  that  pass  through  values  defined  on  the  Gauss- 
Labatto  points.  Since  the  Gauss  points  fall  strictly  between  the  Lobatto  points  [3],  the 
result  is  a  staggering  of  the  solution  and  flux  values. 

The  staggered  grid  leads  to  a  simpler  and  more  flexible  multidomain  method  than 
one  based  on  a  Lobatto  grid  alone,  e.g.  [4].  It  is  simpler  because  subdomain  comers  are 
not  included  in  the  approximation,  so  special  conditions  do  not  need  to  be  derived  for  them. 
It’s  flexibility  comes  from  the  fact  that  only  the  normal  fluxes,  not  the  flux  derivatives, 
need  to  be  continuous  across  interfaces  where  subdomains  meet,  making  grid  generation 
less  restrictive.  See  [1]  for  details. 

The  flexibility  of  the  method  presented  in  [1]  is  still  limited,  however,  by  the 
restriction  that  the  cdculation  of  the  unique  flux  along  a  subdomain  interface  requires  the 
grid  points  to  coincide  there.  We  refer  to  this  approximation  as  conforming.  In  general,  the 
requirement  that  the  interfaces  be  conforming  means  that  subdomains  must  intersect  along 
an  entire  side  or  at  a  comer  point.  If  two  subdomains  intersect  along  a  side,  then 
polynomial  approximation  orders  must  be  the  same  along  the  interface  between  them  [5]. 
An  example  of  a  two  subdomain  conforming  subdivision  of  a  square  is  shown  in  Fig.  la. 

The  limits  imposed  by  the  conforming  restriction  make  it  impossible  to  do  local 
refinement  by  subdividing  existing  subdomains,  or  by  increasing  the  polynomial  order 
within  selected  subdomains.  If  refinement  is  necessary  within  one  subdomain,  it  is 
necessary  to  refine  also  its  neighbors.  This  makes  the  approximation  more  expensive  than 
necessary,  since  the  overall  grid  is  often  refined  where  refoement  is  not  needed. 

To  be  completely  flexible,  we  would  like  the  method  to  be  able  to  use  an  arbitrary 
tiling  of  a  domain  by  quadrilaterals.  This  would  be  similar  to  zonal  finite  difference 
methods  that  have  long  teen  in  use  in  the  finite  difference  community,  e.g.  [6].  Since 
within  each  subdomain  the  strong  tensor  product  structure  of  the  spectral  approximation 
would  remain,  the  result  would  be  a  semi-stractured  method.  The  flexibility  of  the  semi- 
structured  method  would  allow  commonly  available  block  structured  grid  generation 
methods  [7]  to  be  used  to  generate  the  subdomains. 

The  semi-structured  method  can  be  developed  by  loosening  the  restriction  that  the 
fluxes  be  continuous  across  an  interface,  giving  a  non-conforming  patching  of  the 
subdomains.  Fig.  Ib-d  shows  three  non-conforming  topologies.  In  the  first  (Fig.  lb), 
which  we  call  order  refinement,  the  subdomains  intersect  along  a  full  side,  but  the 
approximation  order  changes  across  the  interface.  The  second.  Fig.  Ic,  shows  a  situation 
that  comes  from  a  subdivision  of  the  conforming  grid.  Fig.  la.  Finally,  Fig.  Id  shows  the 
fully  non-conforming  case,  where  the  interface  between  two  subdomains  is  not  a  full  side 
of  either. 
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Fig.  1  Conforming  vs.  non-conforming  grids 

Non-conforming  spectral  domain  decomposition  approximations  for  elliptic 
problems  and  for  the  incompressible  Navier-Stokes  equations  on  grids  like  those  shown  in 
Fig.  1  first  appeared  in  the  late  1980’s  [8-11].  Most  notable  was  the  mortar  element 
method,  in  which  a  one-dimensional  polynomial  function  space  called  a  mortar  was  defined 
along  subdomain  interfaces.  It  was  with  this  mortar  space  that  the  patching  of  the 
subdomains  was  accomplished.  Details  can  be  found  in  the  cited  references. 

In  this  paper,  we  present  a  semi-stmctured  method  that  uses  a  non-conforming 
mortar  approximation  for  the  solution  of  hyperbolic  systems  such  as  the  Euler  gas- 
dynamics  equations.  Interior  to  the  subdomains,  the  method  uses  the  conservative 
staggered-grid  approximation  presented  in  [1].  The  result  is  a  fully  flexible  approximation: 
Subdomains  can  be  subdivided.  Polynomial  orders  can  be  adjusted  within  subdomains 
without  affecting  neighboring  subdomains.  There  is  no  restriction  on  how  the  subdomains 
tile  the  full  domain,  as  long  as  they  do  not  overlap. 

The  paper  begins  with  a  presentation  of  the  equations,  followed  by  a  review  of  the 
staggered-grid  conforming  approximation.  In  Section  4,  we  introduce  the  mortar  method 
for  treating  the  subdomain  interfaces.  Test  problems,  both  linear  and  non-linear,  are 
presented  in  Section  5  to  show  that  the  approximation  is  spectrally  convergent. 
Conclusions  are  drawn  in  the  last  section. 

2.  The  Equations 

In  this  paper  we  consider  the  approximation  of  hyperbolic  systems  in  conservative 

form, 


aQ  aF  ac 

dt  dx  dy 


=  CH, 


(1) 


where  Q  is  the  vector  of  solution  unknowns  and  F(Q)  and  G(Q)  are  the  adveaive  flux 
vectors.  For  the  Euler  gas-dynamics  equations  in  two  space  dimensions. 
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In  (2),  pe  =  p/iy -\)  + p(u^ +  v^)l2  and  we  assume  y  =  1.4.  For  two-dimensional 
problems,  ^  =  0.  For  axisymmetric  problems,  we  interpret  x  as  the  axial  coordinate  and  y 
as  the  radial  coordinate,  and  we  set  1.  For  the  gas-dynamics  equations. 
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In  the  multidomain  approximation  [4],  the  region  under  consideration  is  divided 

into  K  non-overlapping  subdomains,  €1'^.  Each  subdomain  is  mapped  individually  onto  the 
unit  square.  Under  the  mapping,  eq.  (1)  becomes 


dQ  ^  dF(Q)  ^  BGjQ) 
dt  dx  dr 


(4a) 


where  Q  =  7Q  and 


F  =  y;:'F-x^'G  G  =  -y^F-t-x^G  H  =  7H 
J(X,y)  =  x;^yy^-x^y^ 


(4b) 


3.  The  Conforming  Staggered-Grid  Approximation 


The  staggered-grid  approximation  [1]  computes  the  solution  values,  Q,  and  the 
fluxes  F  and  G  on  separate  grids.  These  grids  are  tensor  products  of  the  Lobatto  grid,  Xj, 
and  the  Gauss  grid,  mapped  onto  [0,1] 
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y  =  0,l,...,iV-l 


On  the  Lobatto  and  Gauss  grids,  we  define  two  Lagrange  interpolating  polynomials 
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and 
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We  see  that  £j{x)  eP^yC^)’  and  eP;v-P  where  P^y  is  the  space  of  polynomials  of 

degree  less  than  or  equal  to  N . 

The  mapping  of  each  subdomain  onto  the  unit  square  is  done  by  a  static 
isoparametric  transformation.  Let  the  vector  function  g(s),  0  <  s  <  I  define  a  parametric 
curve.  The  polynomial  of  degree  N  that  interpolates  g  at  the  Lobatto  points  is 

r(5)  =  £g(5y>y(5).  (8) 

J=0 


Four  such  polynomial  curves,  Tmis),  m  =  1, 2,3,4,  counted  counter-clockwise,  bound  each 
subdomain.  As  in  [1],  we  map  each  subdomain  onto  the  unit  square  by  the  linear  blending 
formula 


X^ixj)  =  (1  -  Y)r^(x) + Yr,(X) + (1  -  X)r4(y) + xr^(Y) 

-x,(l  -  X)(l  -Y)-  X^xa  -Y)-  x,XY  -  X4(l  -  X)r  ’ 


(9) 


where  the  xfs  represent  the  locations  of  the  comers  of  the  subdomain,  counted  counter¬ 
clockwise. 

The  solution  unknowns  are  approximated  at  (X,4.,/2,yy+i/2)»z>y  =  0,l,...,A^-l, 
which  we  will  call  the  Gauss/Gauss  points.  The  interpolant  through  these  unknowns  is  a 
polynomial  in  P/y.i.^.i  =  Pat-i  ®  P^-i* 


\-lN-l 

Q(X,K)  =  (10) 

i^O  j^O 

The  horizontal  fluxes  are  approximated  at  the  Lobatto/Gauss  points 
i  =  0,1,..., A;  j  =  0,l,...,iV  - 1 ,  computed  from  the  polynomial  (10) 

)G(e(  ^p^y+1/2  )),  (11) 

Finally,  the  vertical  fluxes  are  approximated  at  the  Gauss/Lobatto  points 

(^i+i/2’lj)’  ^  =  0,1,...,^-!;  j  =  0,1,..., and  are  computed  as 

+1/2’!^;  )G(!2(X, •.1/2.1^,))  (12) 

The  heart  of  the  multidomain  approximation  is  how  the  interfaces  between 
subdomains  are  treated.  In  the  conforming  approximation  (Fig.  la)  the  interface  points 
between  two  neighboring  subdomains  coincide.  However,  the  two  solutions  at  the 
interface  need  not,  since  they  are  computed  independently  from  the  interpolant  through  the 
Gauss/Gauss  points  in  each  subdomain.  From  these  two  values,  however,  a  unique  flux  is 
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computed  so  that  the  characteristic  propagation  of  waves  is  accounted  for.  For  linear 
problems,  we  use  a  flux  vector  splitting.  For  the  Euler  gas-dynamics  equations,  we  use  a 
Roe  solver  [9]  with  the  entropy  fix  to  compute  the  normal  flux  at  the  interface  from  these 
two  values.  Inflow,  outflow  and  wall  boundaries  can  be  treated  by  specifying  the  boundary 
conditions  as  the  extra  solution  value  in  the  Riemann  solver.  Details  can  be  found  in  [1]. 

Once  the  fluxes  are  computed,  we  form  the  semi-discrete  approximation  for  the 
solution  on  the  Gauss/Gauss  grid.  For  each  subdomain 


dt 


+ 

M/2J+1/2 

ax  37 

Ji+l  /  2,j*l  /  2 


li+l/2J+l/2’ 


i  =  0,1,...,A^-1 

7  =  0,1,...,  AT- r 


(13) 


where  the  derivatives ,  defined  as 
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(14) 


are  computed  by  matrix  multiplication.  Equation  (13)  is  then  integrated  in  time  by  a  two- 
level  low-storage  Runge-Kutta  scheme. 


4.  A  Non-Conforming  Mortar  Approximation 

The  only  difference  between  the  conforming  and  the  non-conforming 
approximation  is  how  the  fluxes  are  computed  along  the  interfaces  between  subdomains. 
In  the  conforming  case  (Fig.  la),  there  are  two  solution  values  at  each  interface  point  from 
which  a  single  flux  can  be  computed  directly.  In  the  non-conforming  cases  (Fig.  Ic-d), 
however,  the  grid  lines  do  not  necessarily  match  along  the  interface,  so  that  point  by  point 
transfer  of  information  cannot  be  made  from  a  subdomain  to  its  neighbor. 

We  have  chosen  to  implement  the  transfer  of  information  between  subdomains  by  a 
mortar  method  [9].  The  basic  idea  is  that  the  mortar  (the  “cement”)  cormects  neighboring 
subdomains  (the  “bricks”).  In  our  method,  the  two-dimensional  subdomains  communicate 
only  with  an  intermediate  one-dimensional  construct,  called  a  mortar,  not  with  neighboring 
subdomains  (Fig.  2).  hi  practice,  a  projection  of  the  solution  values  is  made  from  the 
contributing  subdomain  faces  onto  a  mortar.  It  is  on  the  mortar,  and  not  on  the 
subdomains  themselves,  that  the  Riemann  problem  is  solved  to  give  a  unique  flux.  The 
computed  flux  is  then  projected  back  onto  the  subdomain  faces. 

The  use  of  a  mortar  has  several  advantages  over  direct  subdomain  to  subdomain 
communication  of  solution  values.  First,  each  mortar  will  communicate  with  at  most  two 
subdomain  faces.  The  flux  computations  on  a  mortar  can  be  made  independently  of  the 
subdomains  that  contribute  to  it.  Finally,  the  work  of  computing  the  interface  fluxes  is  not 
duplicated. 

The  conditions  for  determining  the  solutions  of  the  hyperbolic  system  along  the 
mortar  are  different  from  those  in  the  elliptic  case  [11].  Since  we  are  computing  the 
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solution  to  the  system  of  equations  in  conservation  form,  we  set  two  necessary  conditions 
to  be  satisfied  by  the  treatment  of  the  interfaces.  The  first  requirement  is  that  the 
approximation  retains  global  conservation.  This  will  determme  the  choice  of  the  projection 
operator  from  the  mortar  back  onto  the  subdomain  faces.  The  second  requirement  we  will 
term  the  outflow  condition. 


Fig.  2  Diagram  of  mortar  communication  between  three  subdomains  that 

subdivide  a  square 

The  outflow  condition  arises  from  the  fact  that,  in  a  hyperbolic  problem,  waves 
should  pass  through  an  interface  unaffected  by  downwind  contributions.  If  the  problem  is 
scalar,  for  instance,  this  means  that  of  the  two  solutions  at  an  interface,  the  solution  from 
the  subdomain  from  which  the  characteristic  comes  (“upwind  side”)  is  used.  To  affect  this 
choice  with  a  mortar,  it  is  necessary  that  the  solution  along  the  upwind  face  be  unchanged 
after  projecting  onto  the  mortar  and  then  back  onto  the  face.  If  this  is  true,  we  say  that  the 
approximation  satisfies  the  outflow  condition. 

We  wUl  describe  mortar  approximations  for  the  three  non-conforming  topologies 
shown  in  Fig.  1.  The  first  two  topologies  occur  when  a  conforming  topology  is  locally 
refined.  The  first  means  only  that  the  polynomial  order  along  a  subdomain  face  differe 
from  that  of  its  neighbor.  The  second  situation  arises  when  a  subdomain  itself  is 
subdivided  without  subdividing  the  neighbor.  The  last  topology  is  the  most  general  one, 
and  does  not  come  from  an  initially  conforming  grid. 

In  the  discussion  that  follows,  we  will  consider  only  the  approximation  of  a  scalar 
problem 


M  +  ^  +  ^  =  0  (15) 

dt  dx  ay 

The  extension  to  a  system  is  direct.  The  only  difference  in  the  case  of  a  system  is  that  along 
the  mortar  a  characteristic  resolution  of  the  two  solutions  must  be  made  to  compute  the 
flux.  That  mortar  flux  calculation  is  identical  to  the  conforming  case  [1]. 
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4.1  Order  Refinement 


The  simplest  case  of  a  non-conforming  approximation  (Fig.  lb)  occurs  when  two 
subdomains  meet  along  a  full  side,  but  the  polynomial  orders  of  the  solution  on  either  side 
of  the  interface  are  not  the  same.  A  diagram  of  two  contributing  sides  of  subdomains  “L” 

and  “R”  and  the  mortar  is  shown  in  Fig.  3.  For  this  topology,  the  mortar,  denoted  by  S, 
extends  the  full  length  of  both  sides  of  the  two  subdomains.  Note  that  the  approximation 
needs  to  consider  only  the  two  subdomains  whose  faces  coincide,  since  the  staggered  grid 
approximation  does  not  include  subdomain  comers. 


L-^H  H  R-^H 

P  ^  P 


(a) 


H, 


R 


R 


(b) 


Fig.  3  Schematic  of  order  refinement,  (a)  Subdomain  to  mortar  projections. 

(b)  Mortar  to  subdomain  projections. 

Only  the  solution  along  the  subdomain  faces  must  be  transferred  to  the  mortar.  So 
define  the  solutions  along  the  faces  as  Uj^i,2=Q^0-,Yj^i,2X  j  =  ^X—,M^ -1 
andC/*^i/2  =  2^(0’^+i/2)’ “1-  Th®  polynomials  along  the  faces  that 
interpolate  these  values  are 


M^-1 

£/"(!)= 

;=o 

C/'(D=: 

j=Q 


(16) 


where  ^  g[0,1]  is  the  local  subdomain  coordinate.  On  the  mortar  itself,  we  represent  the 
two  solutions  as  the  polynomials  (ff  and  defined  by 


(17) 

y=o 

To  be  able  to  satisfy  the  outflow  condition,  the  polynomial  order  of  the  mortar  space  must 
be  sufficiently  large  to  include  both  and  So  as  not  to  require  an  excessive 

amount  of  work  to  compute  the  projections,  we  choose  J  =  max(M^,M*) . 
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To  compute  the  flux  for  each  subdomain  face,  a  three  step  procedure  is  used.  First, 
the  two  solutions  are  projected  onto  the  mortar  space  by  projections  and 

to  give  The  two  values  on  the  mortar  are  then  used  to  compute  a  unique  mortar 
flux,  which  is  evaluated  as  if  the  mortar  is  a  conforming  interface.  Finally,  the  mortar  flux 

is  projected  back  onto  the  subdomain  faces  by  the  projections  P^"^^  and  These 
projection  operators  are  chosen  so  that  the  approximation  is  globally  conservative  and 
satisfies  the  outflow  condition. 


4.1.1  Subdomain  ^  Mortar  Projection 

A  diagram  of  the  projection  of  a  solution  onto  a  mortar  for  the  order  refinement  case 
is  shown  in  Fig  3a.  Since  we  have  chosen  the  order  of  the  polynomial  approximation  on 
the  mortar  (the  “mortar  order”)  to  be  equal  to  the  maximum  order  of  the  polynomials  on  the 
two  contributing  subdomains,  one  of  the  projection  operators  is  tiie  identity.  For 

convenience,  we  will  assume  that  J  =  I^  so  that  and  P®"^®  =  I . 

For  the  projection  of  the  lower  order  space  onto  the  mortar,  we  use  the  unweighted 
projection.  T^us,  we  ask  that 

=  0  m  =  0,1,...,7  - 1  (18) 

Substitution  of  the  definitions  (16)  and  (17)  into  (18)  gives 


^j+l  1 2^m+l  /  2^^  ~  y^^i+l/2[jQ^;+l/2^m+l/2*^^l 

■*  j=0  '■  •' 


m  =  0,1,...,7  —  1 


(19) 


Now  define  the  matrix  elements 


^mj  ~  j^^j+l/2^m+l/2^ 


(20) 


so  that  (19)  becomes  =  MO,  where  U  is  the  vector  of  discrete  solution  values  along 
the  face,  and  O  is  the  vector  of  solution  values  along  the  mortar.  The  integrals  in  (20)  can 
be  computed  exactly  by  a  Clenshaw-Curtis  quadrature  [3]  on  27+2  Lobatto  points.  We 
then  define  the  projection  operator  by 


$  =  =  M-'S^U"-  (21) 

Since  the  computation  of  the  projection  operator  requires  the  inversion  of  the  matrix 
M,  it  is  important  to  consider  the  conditioning  of  that  system.  We  find  numerically  that  the 
growth  of  the  condition  number  in  the  maximum  norm  is  weak  with  matrix  size, 

K  « 1.367°-^.  For  the  maximal  polynomial  orders  that  we  typically  use  (7  =  20),  k  « 15. 
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4.1.2  Mortar  Subdomain  Projection 

Once  the  mortar  functions  (17)  are  computed,  the  normal  flux  is  calculated  as 
described  in  [1].  This  flux  must  then  be  projected  back  onto  the  subdomains  (Fig.  3b). 

Let  ¥  eP,_i  be  the  mortar  contravariant  flux  and  let  F^and  be  the  two  subdomain 

fluxes  to  be  computed  from  it.  Since  =  J,  we  have  immediately  that  =  ¥.  To  get 
the  flux  on  the  left,  we  require  that 

£(f" =  m  =  1,2,.,„M"  -1  (22) 

i.e.,  that  the  projection  error  is  orthogonal  to  ■  Then,  in  matrix  form, 

(23) 

where 


^m,i  ~ 

^m,j  ~  \j^j+\l2^m+\l2^^ 


m,j  =  0,1,...,M^  -1 

m  =  0,l,...,M^  -1 
7  =  0,1,...,/-!  ’ 


As  before,  these  integrals  can  be  computed  exactly  by  quadrature.  From  (23)  the  projection 
of  the  flux  onto  the  subdomain  face  is 


fL  ^  ps^L^  =  (M^  )■*  S¥  (24) 

4.2  Subdomain  Refinement 

The  next  level  of  flexibility  allows  subdomains  be  subdivided  locally.  For 
simplicity,  we  consider  only  the  case  when  a  side  is  subdivided  into  two.  The  refinement 
to  three  or  more  subdomains  is  a  simple  extension  of  this  refinement. 

When  a  subdomain  is  refined  as  shown  in  Fig.  Ic,  there  are  two  possible  choices 
for  the  mortars  (Fig.  4).  In  the  first,  two  mortars  coincide  with  the  “short”  faces  of  and 
The  second  choice  uses  a  single  mortar  that  coincides  with  the  “long”  face  of 


Fig.  4.  Two  mortar  configurations  for  subdomain  refinement 
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We  choose  the  topology  in  Fig.  4a,  since  it  is  the  one  that  can  satisfy  the  outflow 

condition.  In  the  second  case,  the  projection  from  and  onto  the  mortar  is  a 
projection  of  a  piecewise  polynomial  space  onto  a  single  polynomial  space.  The  former  is 
the  larger  space,  since  it  includes  approximations  that  are  discontinuous  at  the  point  where 

and  meet.  The  outflow  condition  requires  that  the  projection  of  face  values  from 

and  onto  a  mortar,  and  the  subsequent  projection  back  onto  the  faces  returns  the  original 
polynomial  functions.  This  is  clearly  impossible  in  case  (4b),  since  the  projection  back 
onto  the  faces  returns  a  continuous  function.  By  using  two  mortars,  however,  as  shown  in 
Fig.  4a,  it  is  possible  to  construct  projections  that  recover  the  original  polynomials  on  all 
three  subdomains.  This  situation  differs  from  the  mortar  element  method  for  elliptic 
problems  [11],  which  require  stronger  regularity  conditions  than  hyperbolic  problems. 

As  before,  we  define  the  solution  approximations  along  a  face  as 


M  -1 


;=0 

M^-1 


(25) 


7=0 

M^-i 


u\S)=  Su'., 2(1) 

j=o 


for  local  subdomain  coordinates  ^  e  [0,1].  We  also  define  four  mortar  functions 


r-1 


7=0 

0"“(z)=iVj^’A;„2(z)sP2.., 


(26) 


j=0 


which  are  functions  of  the  local  mortar  coordinate,  z  €[0,1].  Finally,  we  define  the 
variables  o*  and  /'  to  be  the  offset  and  the  scale  of  a  mortar  with  respect  to  the  subdomain 

U"  that  contributes  to  it.  Thus,  for  z  e  [0,1],  ^  -  o'"  +  s'" z. 

The  orders  of  the  mortar  polynomials  must  be  chosen  sufficiently  high  so  that  the 
outflow  condition  can  be  satisfied.  This  means  that  the  mortar  order  must  be  at  least  as 
large  as  the  largest  subdomain  order  of  all  contributing  subdomains.  Thus,  we  choose 

/'  =  and 


4.2.1  Subdomain  ^  Mortar  projections 

To  compute  the  mortar  functions  (26),  we  also  use  the  unweighted  projection. 
For  each  mortar  S  and  each  subdomain  contributor  Q,  we  require  that 

1 

J (0(z)  -U(o  +  iz))A“.i y  =  0,1,...,/ .  (27) 

0 
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Then  the  vector  of  the  solution  values  along  the  mortar  can  be  computed  by 

$  =  =  M-‘SU  (28) 

where,  this  time 


= 

(29a) 

(29b) 

Note  that  the  matrices  in  (20)  are  just  special  cases  of  those  in  (29)  with  o  =  0  and  5=1. 


4.2.2  Mortar  Subdomain  projections 

Once  the  flux  is  computed  on  the  mortars,  as  if  the  approximation  is  conforming,  it 

is  projected  back  onto  the  subdomain  faces.  The  projections  and  are  exactly 

as  described  by  (22)-(24).  If  the  mortar  order  and  the  subdomain  polynomial  order  are  the 
same,  the  projection  simplifies  to  the  identity,  and  a  simple  copy  of  the  flux  from  the  mortar 
to  the  face  can  be  made. 


The  projection  from  the  mortars  to  a  subdomain  is  a  little  more  complicated  in  the 

case  where  two  mortars  contribute  to  a  subdomain,  as  is  the  case  for  in  Fig.  4a.  In  that 
case,  we  seek  the  flux  that  satisfies 

1 

/  (f(|)  -  (D^l,  j  =  1,2, ...,M  - 1  (30) 

0 


where 


If  we  now  define 


-< 

^  ^  ) 

(i 

^  J 

o'  <^<1 


0 

1 

0 

we  can  write 


(31) 


(32a) 

(32b) 
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(33) 


2 

pi  _ 

*=1 

where  is  the  projection  matrix. 

4.3  Fully  Non-Conforming  Interfaces 

The  final  topology,  shown  in  Fig.  Id,  is  characterized  by  subdomains  whose  faces 
only  partially  overlap.  In  order  to  satisfy  the  outflow  condition,  we  choose  the  mortars  to 
cover  the  intersections  of  two  subdomain  faces,  as  shown  in  Fig.  5.  We  note  again  that 
this  choice  is  different  than  the  choice  one  must  make  in  the  elliptic  case  [9],  where  the 
concept  of  a  supermortar  was  introduced  to  give  the  solution  sufficient  regularity. 


Fig.  5.  Mortar  topology  for  a  fully  non-conforming  approximation 

In  practice,  this  situation  is  handled  exactly  as  in  subdivision.  The  mortar  orders 
are  chosen  to  be  the  maximum  of  the  orders  of  the  polynomials  of  the  contributing 
subdomains.  The  computation  of  the  mortar  functions  is  done  by  eq.  (28)  for  all  the 
subdomains,  since  none  of  the  projection  operators  is  the  identity.  Once  the  mortar  flux  is 
computed  along  each  mortar,  it  is  projected  back  onto  the  subdomain  faces  by  (33),  where 
the  upper  limit  on  the  sum  is  equal  to  the  number  of  mortars  that  contribute  to  the 
subdomain  face. 

4.4  Properties  of  the  Order  Refinement  Projections 

The  use  of  the  unweighted  projections  gives  the  mortar  approximation  two 
desired  properties:  The  approximation  is  globally  conservative  and  the  outflow  condition  is 
satisfied. 

It  is  simplest  to  show  conservation  for  order  refinement.  In  that  case,  we  need  only 
to  show  that 


(34) 

This  is  because  with  ^  =  0  the  integration  of  eq.  (13)  over  all  subdomains  leaves  only 
integrals  of  the  flux  over  the  boundaries.  Conservation  follows  if  the  interface  flux 
contributions  cancel.  By  design. 
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(f‘->I'.v)  =  0  VveP„._, 


and 


(35a) 


(f"'  -  W,w)  =  0  Wwe  (35b) 

where  {u,v)  =  j\vd^.  Since  1  <=  =  Py_„  the  result  follows.  A  similar 

argument  can  be  constructed  for  subdomain  refinement  and  fully  non-conforming 
interfaces.  We  remark  that  (34)  would  not  hold  if  the  mortar  fluxes  are  merely  interpolated 
to  the  subdomain  faces  from  the  mortar. 

The  outflow  condition  is  also  satisfied  by  our  choice  of  projections.  Again,  it  is 
simplest  to  show  this  for  order  refinement.  terms  of  the  projection  operators,  the 

outflow  condition  means  that  =1.  if  we  call  O  =  and  U*  = 

this  is  true  if  U  =  U*.  By  construction, 

(f/-O,v)  =  0Vv  eP^_i  (36a) 


and 


(U  *  -O,  w)  =  0  Vw  e  Pj^i  j  c  Pj_i  (36b) 

Then  (U  -^,w)  =  0  Vw  g  ^  ,  from  which  we  see  that  {U  -  U*,w)  =  0  Vw  g  P^t_i • 
The  result  follows  from  the  fact  that  U  -U*  g  • 

For  the  subdomain  refinement  case,  satisfaction  of  the  outflow  condition  requires 
that  the  projection  operators  satisfy  the  relations 


(37) 

which  can  be  shown  by  arguments  similar  to  that  above. 

4.5  Mortar  Algorithm 

The  algorithm  for  the  non-conforming  approximation  is  the  same  as  the  conforming 
one,  except  for  the  manner  in  which  the  interface  fluxes  are  computed.  At  the  start  of  a 
calculation,  after  the  grid  connections,  mortar  positions,  offsets  and  scales  are  computed, 
the  projection  matrices  are  computed  and  stored.  Then  at  each  stage  of  the  time  integration, 
we  use  the  following  algorithm  based  on  the  staggered  grid  method  of  [1]: 

Algorithm  1.  (Non -Conforming  Staggered  Grid) 

1 . Interpolate  the  Gauss/Gauss  point  solution  values 
to  the  Gauss/Lobatto  and  the  Lobatto/Gauss 
points . 

2. Compute  the  interior  point  fluxes  F  and  G  from 
the  interpolated  values . 
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3 .  Compute  the  Interface  Fluxes : 

(a)  Project  the  interface  solution  values  onto  the 
mortars . 

(b)  Compute  the  mortar  fluxes. 

(c)  Project  the  mortar  fluxes  back  onto  the 
subdomain  faces. 

4.  Compute  the  boundary  fluxes  by  applying  the 
boundary  conditions . 

5 .  Compute  spatial  derivatives  at  the  Gauss/Gauss 
points . 

6 .  Update  the  solution  at  the  Gauss/Gauss  points . 

7 .  Repeat  Steps  1-6  until  done. 

The  mortar  approximation  adds  little  to  the  cost  of  the  calculation  of  the  subdomain 
fluxes  as  N  gets  large,  provided  that  the  projection  matrices  are  computed  and  stored  at  the 
beginning  of  a  computation.  The  bulk  of  the  computations  occur  in  matrix  multiplication 
operations,  which  require  0(A^)  multiplications,  where  N  is  the  order  of  the  matrix. 
Assuming  that  the  same  number  of  points  is  used  in  each  space  dimension,  the  work 
required  to  compute  the  interior  fluxes  for  the  Euler  gas-dynamics  equations  is  0(N^(16N + 
72)).  The  subdomain  faces  are  one  space  dimension  less  and  the  work  required  to  do  step 
(3)  of  the  algorithm  above  is  0(N(16N  +  132)).  Thus,  the  work  required  to  treat  the 
interfaces  relative  to  the  interior  work  is  asymptotically  0(l/N). 

5.  Examples 

In  this  section,  we  use  the  semi-structured  algorithm  to  compute  solutions  to  both 
linear  and  non-linear  hyperbolic  problems.  We  first  solve  a  two  variable  linear  system 
using  non-conforming  topologies  and  compare  the  convergence  to  the  convergence  using 
alternative  conforming  grids.  We  then  present  an  example  where  the  solution  is  localized, 
and  show  that  the  computational  cost  for  the  same  error  can  be  reduced  significantly  by 
using  the  non-conforming  interface  treatment. 

We  also  solve  three  problems  using  the  Euler  gas-dynamics  equations.  The  first 
problem  is  that  of  subsonic  flow  from  a  point  source,  for  which  there  is  an  exact  solution. 
The  convergence  rates  of  conforming  and  alternative  non-conforming  grids  are  compared. 
We  then  show  that  exponential  convergence  is  obtained  when  solving  the  problem  on  a 
complex,  multiply  connected  subdomain  topology.  The  second  problem  is  a  steady 
subsonic  flow  over  a  circular  bump.  We  show  that  exponential  convergence  of  the  entropy 
is  obtained  for  both  conforming  and  non-conforming  approximations.  The  non-conforming 
approximation,  however,  takes  only  half  the  computer  time  for  the  same  accuracy.  Finally, 
as  an  example  of  a  transonic  flow  problem  we  solve  flow  in  an  axisymmetric  converging- 
diverging  nozzle.  The  computed  results  for  that  problem  are  compared  to  experimental  data. 

5.1  Linear  Model  Problem 

We  begin  by  considering  steady  solutions  to  the  system 

Q,-hF,-f-G,=S(x,3;)  (38a) 


where 
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Q  = 


F  = 


3v-u 
3m- V 


,  G  = 


2m  +  4v 
4m+2v 


The  source  term  is  chosen  so  that  the  exact  steady  solution  is 

U  =  ^sm(2jt;x)sin(2iiy) 


(38b) 


Exponential  convergence  is  observed  for  all  three  non-conforming  grid  topologies 
shown  in  Fig.  lb.  First,  we  consider  order  refinement  on  a  subdivision  of  the  rectangle 

[0,2]  X  [0,1].  Fig.  6  compares  the  errors  of  a  conforming  grid  with  two  subdomains 
with  those  of  a  non-conforming  approximation.  In  both  cases  we  observe  exponential 
convergence.  For  order  refinement  alone,  we  would  expect  the  error  to  be  dominated  by 
the  lowest  order  approximation,  and  this  is  the  case. 

Exponential  convergence  is  also  observed  when  a  subdomain  is  subdivided.  Fig.  7 
compares  the  error  of  a  four  subdomain  conforming  decomposition  with  a  subdivided 
approximation.  As  before,  we  observe  that  the  error  is  dominated  by  die  A""  order 
polynomial  approximation. 

Finally,  we  consider  a  fully  non-conforming  subdivision,  shown  in  Fig.  8.  In  this 
case,  we  choose  a  subdivision  Aat  is  between  two  conforming  approximations  in  its 
resolution.  Again,  we  observe  exponential  decay  of  the  error,  and  that  error  lies  between 
the  errors  of  the  two  conforming  approximations. 


N 


Fig.  6.  Comparison  of  conforming 

and  non-conforming  errors  for 
order  refinement. 


N 


Fig.  7  Comparison  errors  for 

conforming  and  non-conforming 
grids. 
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-1 

-2 


N 


Fig.  8.  Comparison  of  conforming  and 
fully  non-conforming  errors. 


(a)  Solutions 


(b)  Conforming  Grid 


(c)  Non-conforming  Grid 


Fig.  9.  Solution  and  grids  for  (38)  where  local  refinement  is  needed. 

One  of  the  main  reasons  to  use  a  non-conforming  grid  is  to  compute  efficiently 
solutions  where  local  refinement  of  the  grid  is  needed.  As  an  example,  we  compute  the 
solution  of  (38)  on  the  unit  square  where  the  source  terms  are  chosen  so  that  the  steady 
solution  is 


v  =  e 


-50{x^+y^) 


In  this  problem,  the  solution  variations  are  concentrated  in  the  upper  left  and  lower  right 
comers  as  shown  in  Fig.  9a.  We  solve  the  problem  on  two  grids,  also  shown  in  Fig.  9 . 
The  non-conforming  grid,  which  has  increased  resolution  only  where  needed,  has  44%  of 
the  number  of  grid  points  of  the  conforming  grid  when  N  =  10. 

Fig.  10  compares  the  convergence  of  the  error  for  die  two  grids  shown  in  Fig.  9. 
We  find  that  the  errors  for  both  grids  are  the  same  to  one  digit.  However,  at  A  =  10,  the 
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computational  cost  of  the  non-conforming  grid  computation  is  46%  of  the  cost  of  the 
conforming  grid  computation.  Since,  at  best,  we  would  expect  the  non-conforming  grid  to 
take  44%  of  the  time  of  the  conforming  one,  we  see  that  the  oveihead  due  to  the  mortar 
projections  is  only  about  5%. 
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Fig.  10.  Convergence  of  linear  model  problem 
(38)  on  the  two  grids  shown  in  Fig.9. 


5.2  Euler  Gas-Dynamics  Equations 

5.2.1  Point  source  flow.  We  now  consider  the  solution  of  the  flow  of  a 
steady,  iirotational  gas  exiting  from  a  point,  which  can  be  solved  exactly  by  a  hodograph 
transformation  [13].  The  streamlines  are  radial  and  the  level  curves  of  the  Mach  number, 
pressure  and  density  are  circles  centered  on  the  source.  We  wiU  compute  this  flow  in  two 
geometries.  The  first  geometry  represents  flow  in  an  expanding  duct,  where  two 
streamlines  are  chosen  as  walls  of  the  duct.  The  second  geometry  is  that  of  a  rectangular 
region  with  three  cut  out  circles. 

The  first  geometry  represents  steady  flow  in  an  expanding  two-dimensional  duct 
with  straight  walls  (Fig.  11).  The  lower  wall  was  chosen  to  be  the  line  y  =  0  and  the  upper 
wall  the  line  y  =  x  tan(7c/6).  The  exact  solution  chosen  was  the  one  that  takes  on  the  Ji^ch 
number  M  =  0.6  at  the  lower  left  comer.  We  compute  this  flow  on  the  two  grids  shown  in 
Fig.  11.  An  examination  of  the  error  using  a  single  domain  approximation  indicates  that 
most  of  the  contribution  of  the  error  occurs  near  the  lower  left  comer.  Thus,  we  set  up  the 
non-conforming  grid  as  shown  in  Fig.  1  lb.  For  comparison,  we  also  compute  the  solution 
on  the  conforming  grid  shown  in  Fig.  1  Ic.  The  errors  for  the  two  grids  are  plotted  m  Fig. 
12,  and  are  the  same  to  one  digit. 
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(a)  Mach  Contours  (b)  Non-conforming  grid  (c)  Conforming  grid 

Fig  11.  Diffuser  solution  and  non-conforming  and  conforming  grids. 
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Fig.  12  Convergence  of  the  density 
error  for  the  two  grids  of  Fig.  11. 


Exponential  convergence  can  be  obtained  on  complex  geometries,  too,  if  the 
solution  is  smooth.  In  Fig.  13,  we  show  the  Mach  contours  and  grid  for  the  solution  of  the 
point  source  flow  in  a  rectangular  domain  with  three  cut  out  holes.  The  point  source  was 
placed  at  the  center  of  the  bottom  circle.  This  grid  has  fully  non-conforming  subdomain 
interfaces.  The  exact  solution  was  used  to  compute  the  boundary  conditions  on  all 
boundaries.  Fig.  14  shows  that  the  error  converges  exponentially  with  N. 
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Fig  13.  Mach  Contours  and  multiply  connected  grid  for  the  point  source  flow. 


Fig.  14.  Convergence  of  the  density 
for  the  grid  in  Fig.  13. 

'2  4  6  8  10121416 
N 

5-2.2  Subsonic  Flow  Over  a  Circular  Bump.  Our  next  example  is  the 
solution  of  a  Mach  0.3  subsonic  flow  over  a  circular  bump.  This  flow  was  computed  on 
two  grid  topologies,  shown  in  Fig.  15.  A  wall  boundary  condition  was  specified  dong  the 
bottom.  At  the  left,  right  and  top  boundaries,  the  uniform  free-stream  condition  was 
specified  as  the  external  input  to  the  boundary  Riemann  problems.  Initially,  the  free-stream 
solution  was  specified  everywhere,  and  then  the  boundary  conditions  were  imposed.  A 
comparison  of  the  pressure  computed  on  the  conforming  and  non-conforming  grids  is  also 
shown  on  Fig.  15. 
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This  problem  does  not  have  an  exact  solution.  However,  since  the  free-stream  is 
irrotational  and  homentropic,  the  entropy  must  remain  constant  everywhere.  That  this  is 
not  true  computationally  is  due  to  spatial  approximation  errors.  Fig.  16  shows  the 
exponential  convergence  of  the  entropy  error  for  the  two  grids  shown  in  Fig.  15. 


Fig.  15.  Pressure  contours  and  grids  for  M  =  03  flow  over  a  circular  bump.  The  pressure 
contours  for  the  corforming  grid  are  plotted  with  dashed  lines;  the  non- 
conforming  solution  is  plotted  with  solid  lines. 
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N 


Fig.  16.  Convergence  of  the  entropy 
for  the  two  grids  shown  in  Fig.  15. 


For  A^=  10,  the  non-conforming  grid  in  Fig.  15  has  50%  of  the  number  of  degrees 
of  freedom  of  the  conforming  grid.  We  find  that  the  CPU  time  for  the  non-conforming 
approximation  to  get  to  steady-state  is  44%  of  the  time  required  by  the  conforming  one  so 
the  work  required  by  the  mortar  projections  is  negligible. 

4.2.3  Transonic  Flow  in  a  Converging-Diverging  Nozzle.  Our  final 
example  is  that  of  a  transonic  flow  in  an  axisymmetric  converging  diverging  nozzle.  We 
use  the  nozzle  of  Cufifel  et  al.  [14].  The  nozzle  has  a  converging  section  with  half  angle  of 
45°  and  a  diverging  section  of  15°.  The  experimental  tests  were  done  in  air  with  a 
stagnation  temperature  of  540  R  and  a  stagnation  pressure  of  70  psia.  The  nozzle  geometry 
and  grid  are  shown  in  Fig.  17.  We  have  increased  the  resolution  in  the  neighborhood  of 
nozzle  wall  curvature  singularities  by  subdivision  of  a  conforming  grid  (See  [1]). 
Boundary  conditions  and  scaling  were  treated  as  in  ref.  [1]. 
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Fig.  17  Nozzle  shape  and  grid  for  transonic  flow  computation 


Fig.  18.  Computed  and  measured  (symbols)  Mach  contours  in  the  nozzle. 

Results  computed  for  the  nozzle  are  shown  in  Figs.  18-20.  Contours  of  the  Mach  number 
are  compared  to  the  experimentally  determined  positions  in  Fig.  18.  Wall  values  of  the 
pressure  are  and  Mach  number  are  shown  in  Figs.  19  and  20  for  different  subdomain 
resolutions. 
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Fig.  20.  Wall  Mach  number  for  the  converging-diverging  nozzle. 
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6.  Summary 

It  was  difficult  to  do  local  refinement  of  the  grid  with  the  original  staggered-grid 
multidomain  approximation!!].  That  method  required  subdomains  to  intersect  either  along 
a  full  side  or  at  a  single  point.  In  addition,  along  adjoining  faces  the  approximation’s 
polynomial  orders  had  to  be  the  same.  This  made  it  impossible  to  subdivide  a  subdomain 
or  to  increase  locally  the  polynomial  order  as  necessary  to  resolve  a  local  feature  in  the 
solution. 

In  this  paper,  we  have  described  a  semi-structured  method  that  uses  the  staggered 
grid  scheme  interior  to  subdomains.  It  relaxes  the  restriction  that  the  fluxes  be  continuous 
at  subdomain  interfaces  and  allows  for  non-conforming  interfaces.  This  makes  it  possible 
to  use  a  general  quadrilateral  tiling  of  a  computational  domain.  Subdomains  can  be 
subdivided  as  necessary.  Within  each  subdomain  the  polynomial  order  can  be  set 
independently  of  its  neighbors.  More  generally,  subdomains  need  intersect  only  partially 
along  a  side. 

In  the  semi-stmctured  version,  the  interfaces  are  treated  by  a  mortar  method.  The 
solutions  along  subdomain  faces  are  first  projected  onto  a  one  dimensional  construct  called 
a  mortar.  Along  the  mortar  a  unique  flux  is  computed  as  if  the  approximation  is 
conforming.  That  flux  is  then  projected  back  onto  the  subdomain  faces  to  be  used  to 
update  the  solutions  within  the  subdomains.  Asymptotically,  the  work  involved  in  this 
process  is  small  when  compared  to  the  work  required  to  update  the  solutions  in  the  interiors 
of  the  subdomains.  The  mortar  approximation  is  designed  so  that  the  method  is  globally 
conservative. 

The  method  has  been  applied  to  both  linear  and  non-linear  smooth  problems  for 
which  exact  solutions  are  known.  In  all  cases,  spectral  convergence  was  observed.  The 
advantage  of  being  able  to  locally  refine  the  grid  was  also  seen  in  the  reduced  cost  of  the 
non-conforming  approximations. 
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